function cubic_spline_interpolation()
    % 定义函数f
    f = @(x) 1./(1 + 25*x.^2);
    
    % 插值点数N的数组
    N_values = [6, 11, 21, 41, 81];
    errors = zeros(length(N_values), 1); % 存储最大误差
    
    % 准备绘图
    figure;
    hold on;
    colors = ['k', 'r', 'b', 'g', 'm']; % 颜色数组，用于不同的N值
    
    % 绘制真实函数
    xx = linspace(-1, 1, 1000);
    plot(xx, f(xx), 'k-', 'LineWidth', 2);
    
    % 对于每个N值进行操作
    for i = 1:length(N_values)
        N = N_values(i);
        
        % 生成节点和计算真实函数值
        x = linspace(-1, 1, N);
        y = f(x);
        
        % 创建立方样条插值
        pp = spline(x, [0 y 0]);
        
        % 计算中点和误差
        x_mid = (x(1:end-1) + x(2:end))/2;
        y_spline_mid = ppval(pp, x_mid);
        y_exact_mid = f(x_mid);
        errors(i) = max(abs(y_spline_mid - y_exact_mid));
        
        % 在同一图形上绘制插值样条
        plot(xx, ppval(pp, xx), '--', 'Color', colors(i), 'LineWidth', 1.5);
    end
    
    % 完成绘图设置
    title('Cubic Spline Interpolation for Different N');
    legend(['Exact Function', arrayfun(@(N) ['N = ' num2str(N)], N_values, 'UniformOutput', false)]);
    xlabel('x');
    ylabel('f(x)');
    hold off;
    
    % 输出误差和收敛率
    disp('N values:');
    disp(N_values');
    disp('Maxnorm of interpolation errors:');
    disp(errors);
    if length(errors) > 1
        rates = log2(errors(1:end-1)./errors(2:end));
        disp('Convergence rates:');
        disp(rates);
    end
end

